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A Divide and Conquer Method for Unitary and Orthogonal E i genprob 1 ems 
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L . Re i che 1 
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Abstract: Let H G be a unitary upper Hessenberg matrix whose eigen- 

values, and possibly also eigenvectors, are to be determined. We describe 
how this eigenproblem can be solved by a divide and conquer method, in 
whicli tile matrix H is split into two srnallei' unitary right Hessenberg 
matrices II and W 2 by a rank -one modification of H. I'he e i genp rob 1 ems for 
and II 2 can be solv^ed independently, and the solutions of these smaller 
e i genprob 1 ems define a rational function, whose zei^os on the unit circle 
are the eigenvalues of H. I'he eigenvector's of II can be determined from 
the eigenvalues of II and the eigenvectors of llj and II 2 . The outlined 
splitting of unitary upper Hessenberg matrices into smaller sucli matrices 
is carried out recursively. Ihis gives r^ i se to a divide and conquer 
method that is suitable for implementation on a parallel computer^. 

When 11 E is ortliogonal , the divide and conquer scheme simp] if ies 

and is described separatel^^. Our^ inter^est in the ortliogonal eigenproblem 
stems from applications in signal processing. Numerical examples for the 
orthogonal eigenproblem conclude the paper. 

Subject classification: AMS(MOS): 65F15: CR : 1.3 

Key Words: unitary eigenproblem, orthogonal eigenproblem, divide and 

conquer, parallel algorithm, Pisarenko frequencies, Gauss-Szego 
quad ratu re . 
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1 . I nt^roduct ion 

Divide and conquer (DC) methods have been developed for the symmetric 
e i genprob 1 em , see Cuppen [Cu] , Dongarra and Sorenson [DS] , and 
Krishnakumar and Mort [KM] , and have for these problems been shown to be 
efficient on parallel computers and competitive on single processor 
machines [DS] , [Cu] , [KM] . The DC method has also been applied successfully 
to the computation of singular values by Jessup and Sorensen [JS] . In the 
present paper we describe a DC method for the unitary e i genprob 1 em , and we 
also discuss the simplifications that arise for real orthogonal matrices. 

Let H E be unitary. Then H is unitarily similar to a upper 

Hessenberg matrix with real-valued non-negative subdiagonal elements. If 
a subdiagonal element vanishes, then the eigenproblem splits into 
e igenprob 1 ems for smaller upper Hessenberg matrices. We therefore may 
assume that H is a upper Hessenberg matrix with positive subdiagonal 
elements. Then all eigenvalues of H are simple. It is easily seen that H 
can be written as a product of n Givens reflectors Gj E 

H = H(7^,729.*«97n)*=G^G2«*‘G^_jGn9 (1.1) 

where, for 1 < k < n, 






■^k ^k 
^k Tk 

^n-k-1 



7k ^ E IR , > 0 , 

I ')'k I ^ 1 ’ 



and 



Gn : = 



■ n-l 



-7n 



; 7n e C , I I = 1 



(1.2a) 



(1.2b) 
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Here Ij denotes the j x j identity' matrix. The 7 j , 1 < j < n, are the so- 

called Sch u r parameters of H, and 7 j denotes the complex conjugate of 7 j . 

The cr j , 1 < j < n, are said to be compl ernentary parameters of II, and are 

the subdiagonal elements of H. 

The DC method described uses the product representation ( 1 . 1 ) of H, 
the so-called S c h u r parametric form of II . An application of particular 
interest to us is the computation of Pisarenko frequency estimates for a 
random stationary stochastic process, see below. In this application H is 
defined by its Schur parametric form. The determination of Gaussian 
quadrature rules on the unit circle, so-called Gauss-Szego quarature 
rules, also gives rise to unitary (or real orthogonal) matrices in Schur 
parametric form, see Section 5. When the Schur parameters are not 
explicitly known, they can be computed from 



Tl = 


-IIll 






V = 


.pll pH 
- (Gj_i Gj 2 . 


• • ^2 “)jj’ j = 2,3, . . 


, . , n , 



where denotes the conjugate transpose of G|^ , and Mjj denotes the element 

(J^j) of a matrix M € . 

The DC method is most easily described for II G orthogonal. Then 









•'l 




Is 


H = Gi G2 . 


• • Gj.i Gs Gs + i . 


. . Gn = : 


i n-s 


Gs 


H2 



where II^^ G and II2 G are orthogonal. The Giv^ens i^ef lector 

Gs € R^^^ , s < n, can be written as a Householder transformation 

Gs = I - 2 ww^ (1-4) 

* where : 

I 
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w : = 




(1.5) 


u;s : = 


+ 7s) , 


( 1 . 6a) 




;= - 7s) . 


(1 .6b) 



Throughoul: this paper ej denotes the jth column of an identity matrix of 
appropriate order. By (1.3)-(1.4), H is orthogonally similar to 



h' 




(I - 2w^) 



fi - 2Hww^ . 



(1.7) 



This is one step of the DC method for orthogonal matrices: The eigen- 

values, and if so desired the eigenvectors, of and H 2 are computed 
first. is a rank-one modification of H, and the eigenvalues of H are 

computed as the zeros on the unit circle of a rational function, whose 
poles are the eigenvalues of . 

Section 2 describes the DC method for unitary matrices. In Section 3 
we show some results on the orthogonality of the eigenvectors and on the 
location of the eigenvalues. These results are analogous to bounds 
presented by Dongarra and Sorensen [DS] for the DC method for symmetric 
matrices. Section 4 discusses simplifications that can be made when 
H is real and orthogonal, and also considers some computational details. 
Computed examples for the orthogonal eigenvalue problem are presented in 
Sect ion 5 . 

An outline of a unitary DC method with convergence results for the 
root-finder has previously been presented in [GR] . The splitting into 
subproblems is done differently in the present paper. A related DC method 
is described by Arbenz and Golub [AG] . Cybenko [Cy] reduces the 
orthogonal eigenproblem to an eigenproblem for a symmetric tridiagonal 



4 



matr i x . 



The orthogonal e i genprob 1 ern is in [AGRl] solved by solving 



singular value problems for certain bidiagorial matrices, and a [Jll 
algorithm Toi" unitary matrices is presented in [Grl] . A comparison with 
respect to accuracy and speed of these methods still remains to be done. 
Here we only note that DC methods are suitable for implementation on 
parallel computers, see [DS] , [KM] , and Section 2. Some of the schemes 
mentioned transform the orthogonal eigenvalue problem to a symmetric one. 
The real eigenvalues of the latter problem are then map^ped to the unit 
circle to yield the eigenvalues of the orthogonal e i genprob 1 em . The 
mapping from the interval to the unit circle may be sensitiv'e to 
perturbations of arguments near the end points of the interval, and it may 
therefore be difficult to determine eigenvalues close to ±1 accurately 
with these schemes. 



Pisarenko [Pi] proposed a method for decomposing a random stationary 
stochastic process q ^ ^ into a sum of harmonics in white 

no i se , i . e . , 

P 

Xm = E cos(m?!>« + . tn > 0 , (1-^) 

«=1 



where the 0^ are arbitrai'y pliase shifts and {ym}|^_Q i® zero mean white 
noise process with variance a . The <f)^ are called P i saren ko frequency 
est i mates . Assume for simplicity that p is the number of distinct 
harmonics in the ^signal’ {x, 7 i}^_q is known, and that 0 < (f>^ < n for 
1 ^ ^ ^ P* Tlien the <j>^ can be determined as follows. Form the 

(2p+l ) X (2p-f 1 ) Toeplitz covariance matrix M for the signal {xm)[^_Q . and 
compute its least eigenvalue • Then ^ see [Pi] • Let be 

the Schur parameters associated with the Toeplitz mati ix M-A^l^l. They can 
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be determined from the Szego recursions (Levinson’s algorithm), see e.g. 
[AGR2] . From our assumptions it Follows that is singular, but 

leading principal submatrices not identical with are not. 



Tlierefore, - 1 < 7 j < 1 for 1 < j < 2p , and 72 p G {-1,1}. By (1.1)-(1.2) it 

follows that the Schur parameters { 7 j}^^j define an orthogonal matrix H G 



IR 



2PX2P 



.2p 



with distinct eigenvalues {Aj}j_j . Enumerate the eigenvalues so that 



those with Im(Aj) > 0 have smaller index than the eigenvalues with Im(Aj) < 
0. Then the Pisarenko frequency estimates are given by 



:= arg(Aj) , 1 < j < P • 

The coefficients Oj of (1.8) are two times the weights belonging to the 
Gauss-Szego quadrature rule with abscissas Aj , 1 < j < p* For details see 

[AGR2] , where also references to related work can be found. The unitary 
DC method yields the Gauss-Szego weights with no extra computational 
effort when computing the eigenvalues Aj . Gauss-Szego quadrature is 
discussed in Example 5.1 of Section 5. 
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2 . 



The unitary e i genprob 1 ern 



In this section we describe a divMde and conquer method for unitary 
right Hessenberg matrices with positive subdiagonal elements. First we 
need to generalize the splitting (1.3)~(1.7) to Givens reflectors with 
compl ex-val ued Schur parameters. This is accomp 1 i slied by noting that tlie 
Gs defined by (1.2a) are diagonally unitarily equivalent with a real 
Householder transformation. Introduce 

f 7s/ I 7s I » 7s 7^ 0 

rS := 

1. 1 . 7. = 0 . 



Then | 7s I = 1? and Tor Gs defined by (1.2a) we obtain 



'is+l 




Is 


7s 


Gs 


7s 


In-s 




^ n-s-1 



^ s-1 



-|7s I <^s 
^s I 7s I 



n-S“l 



( 2 . 1 ) 



where, similarly to (1.5)-(1.6), w = CsWg + ^s+l^^s+l ^ 



Ws := + hsl)'^^ 



u-3+1 := -2'/"(l - hsl)'^^ 



(2.2a) 

( 2 . 21 )) 



Substitution of (2.1) into (1.1) yields 



H, 



H 



In-s 



(I - 2ww^) 



(2.3) 



with 



Hi : — H (7i , 72 ) • • • ) 7s-i ? 7s ) ? 

»2 == H(7sSs+1 > TsSs+ 2’ • • • ’"^sSn) • 



A unitary similarity transform of (2.3) yields, analogously with (l.r), 



H' : = 



lli 



H. 



(I - 2ww*^) =: H - 2Hww^ 



(2.4) 



Let 



Hk = Wk Ak , k = 1, 2, 



(2.5) 



be spectral resolutions, i.e., the are unitary and the are diagonal. 

Then H has the spectral resolution H = W A , where 









Ai 






— 


W2 


, A := 




= diag(Aj,A2, . . 


• 5 ^n) 9 



with I I = 1 for 1 < k < n. 

We are in a position to describe how the spectrum of H can be obtained 
from A, the last row of Wi and the first row of W 2 . Introduce the 
characteristic polynomial 

X(A) := det(AI-H) = det(AI-H') = det (AI -fl+2Hww'^) 

= det(Al-H) det(I+2(AI-H)'^Hww^) 

= + 2w*^(AI-fl)'^fiw) 

= t/’(A)( 1 + 2w^W(Al-A)‘^AW^w) , 
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where 11^ is defined by (2.4), W and A by (2.6), w by (2.2) and V’(^) • = 

def(A[-H]). Let z = be given by 



z : zz W*^w 



Wj’eso-s 

^2 



(2.7) 



and define the spectral function 



-A(A) := = 1 + 2z^^(AI-A)-1Az = 1 + 2 £ Kj | ^ 



n « A -|- Aj 

= E Kjl" 

j=i ^ 



( 2 . 8 ) 



|_i 

where we have used that z z = 1 . Let 



:= arg(Aj) , 0 := arg(A) , 0 < ^ j , ^ < 



27T 



Then, with i := ^ - 1 , 



^(A) 

4 >'( 0 ) 



i ,E ICjl^ cot^-i-2 — ) =: i<i>(6») , 



1 

2 * 



(2.9) 



( 2 . 10 ) 



We may assume that the 0^ are distinct and that all (j ^ 0, because 

otherwise we can make these conditions hold by deflation, see below. Let 
^j' e [o,27t[, 1 ^ j denote the zeros of 4>(^?). Then the eigenvalues of 

and of II are given by Aj^ := exp(i6^j^) , 1 < j < n. The sets and 

strictly interl ace . 

We describe a rootfinder for 4>(^). By the inequality (2.10), the 
zeros of 4>(^) can be determined accurately. We may assume that 
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0 < 6 ^ < 62 < • • • < < 27t and that ^ our initial approximation of a 

zero of satisfies ^n- 2 ;r < < 6 ^, By the strict interlacing of the 

sets and has precisely one zero, denoted , in the 

open interval ]0n-27r, 0^\_, Assume for the moment that 



< 0 , 



( 2 . 11 ) 



and introduce 



- /S^ - 6\ 

^{ 6 ) := p a cotf ^^-2 j • 



( 2 . 12 ) 



The coefficients p and a are determined by osculatory interpolation, i,e.. 






(2.13) 



which yields 



'p = $( 0 ^°^) - 4>'(0^°^) sin((9i - 0 ^°^) , 






.cr = 2 <t>'( 0 ^°^) sin 2 (^ 



) • 



The zero of in ]^n“27r,^j^[ is our next approximation of 6 ^^, New 

approximations 6 ^^^ ^ of 0 ^^ are computed from j > I 9 in a similar 

fashion. The sequence satisfies 9^^^ < 9^^ for j >0, and converges 

monoton i cal 1 y and quad rat i cal ly to 9^^ as j increases, see [GR] for a 
proof , 

If instead of (2.11) we have 



$(^^°^) > 0 , 



(2.14) 



then we replace ( 2 . 12 ) by 



^( 0 ) = p + (T cot^^l!- 2 — 



(2.15) 
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in (2.13). This defines p and a. The zero 6^^^ of 4>(^?) in t.he open 
interval ]^n“27r,^j^[ is our next approximation of 0^^ . New approximations 
of 0^^ are computed from for j > 1 in a similar fashion. The 

sequence satisfies 6^^^ > 0^^ for j >0, and converges monoton i cal 1 y 

and quadrat i cal ly to , see [GR] . In the implementation used to generate 
the computed examples of Section 5, the iterations are carried out until 
^(^(j+l)) ^ The value 0^^^ is accepted as an approximate 

root of = 0. 



F rom 



A : = d i ag [< 






\ 6 ^ 



> 0 ,,'. 



(2.16) 



and the spectral resolutions (2.5) of \\^ and H 2 > we can now compute the 
spectral resolution of H: 



H = WAW^ , W'^W = I . 



(2.17) 



By (2.3)-(2.4) , we can Tor some vector u € C" express II as 



H = fi - 2uu'^ , 



H^est 






6 C" 



Let A := exp(i0^) and v 6 C” form an eigenpair of II, i.e. llv = vA . Then 



(II - 2uu*^)v = vA , 



or, equivalently, 

(n “ IA)v = ua , Q : = 2u*^v . 



This shows that any normalized eigenvector v' of H associated with the 
eigenvalue A is a normalization of 



( 2 . 18 ) 



V 



/ 



(Hl-lA)-^Hiesc^s 

(H2-IA)-leiW5+l 

Wi(I-A}'A)'Hv}‘esWs 
W2(A2-IA)-^v”esW3^1 ’ 



where Wj^ and are given by (2.5). Leb || H 2 denote the Euclidean vector 

and matrix norms. From ||W ^||2 = Ih^2ll2 “ ll^ilb ~ ^ and (2.6), (2.7), (2.10), 

(2.18) it Follows that 

6(A) := ||v'||2 = l|(A-IA)-l ^2 

= (E > J • (2-19) 

j=l |Aj - A| ^ 



We choose 



^A 



Wi(I-A5^A)-^w5*esWs 

W2(A2~IA) 



/KA) , 



(2.20) 



and note that the lower bound (2.19) For <5(A) indicates that severe 
cancellation oF signiFicant digits does not take place in the computation 
oF v^ by (2 . 20) . 

By (2.7) , we only require the last row oF VI ^ and the First row oF W 2 
(as well as A) in order to determine the spectral Function (2.8) . Hence, 
iF we do not desire the eigenvectors, then it suFFices to determine the 
First and last rows oF W in order to be able to compute the spectrum oF 

IT U 

larger problems. We call the triplet {A,e"W,enW} the partial spectral 
resolution oF H. The First and last elements oF v^ can easily be 
determined From 



je«v^ = efWi(I-A«A)-lwHesa;s/KA) , 
lej?v^ = eLw2(A2-lA)-lW»eiW3+l/6(A) . 
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We may assume that the columns of W are such that al 1 corn})onents c:> f the 
vector W*^ej^ are real and positive. Then ej^Uv*^ej^ is tiie weig^tit corresj>ond- 
ing to tile node exp(i^j^^) in the Gauss-Szego quadrature rule with nodes 
{exp ( i , see [Gr2] and Example 5.1. 

We assumed above all components Q of z to be non - van i sii i ng and all 
eigenvalues Aj of H to be distinct. These conditions can be made to iiold 
true by de f 1 at ion . Our discussion follows Dongarra and Sorensen [DS] . 

First assume that Q vanishes. By (2.4)-(2.7), 

W^n'w = A - 2W’^fiww'^W = A(I - 2zz^) , (2.22) 

and from z = 0 it follows that 

A(I - 2zz*^)e£ = Ae^ = * (2.23) 

Substituting (2.23) into (2.22) yields 

H'We^ = We^A^ , 
and therefore 



1 

t 










- 




WjAj^ 



Thus if = On then we can determine an eigenpair of 11 without oxf>licitl\' 
computing a zero of ^(0) and without using (2.20). 

For a general z 6 , with z^z = 1 , we obtain 

||A(I-2zz^)eg||2 = 2U^\ , 

and we a,ccep"t {Ag,e^} as an eigenpaii' of A if 
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(2,24) 



2|C£| < Cl 

f^or some smal 1 constant: , 

Assume that z = [Cj]jLi with 2 | Cj I > al 1 j , We may now be able to 

deflate due to close eigenvalues. Let A = diag[A^,A2, . . . ^ An] with 
Al « A2 , and choose the Givens reflector 

G , (2.25) 

so that 

Gz = [(Kil^ + IC2l^)^^^0^C3»C4 >Cn]'^ , 

i . e . 

:= Ksl/CKil" + , 

7 ••= -Cl |-^/(ICll" + l<2l")'^" • 

\v'e accept + 72 I 7 1 ^ > G*^G2} as an (approximate) eigenpair of 

A(I - 2 zz^) if 

|t<t(Ai - A2) I < ^2- ( 2 . 26 ) 

for some smal 1 constant e2 5 because 

||A(I- 2 zz^)G'^e 2 - (Ai(t2 + 72 I 7 1 ^)G^e 2 ll 2 = l 7 <^( 7 i - 72) I • 

If = 72 then we have determined an eigenpair exactly. In case 
A^ A2 we note that |7cr(A^ - -^2) I — ^ I " -^2 I ’ and, moreover, if I7I « 0 

or I7I « 1 then |7cr(A^ - A2) | << \ - A2 | . Hence, inequality ( 2 . 26 ) may 

be satisfied, even if \ - A2 | > ^2 • Assume that ( 2 . 26 ) is valid. Then A 

is replaced by 
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(n-i)x(ii-l) 



A : — d i 3-g C "T ) ^2^ 9 ^3 ) 9 • • * 9 '^nl ^ ^ 

and if A has close eigenvalues, then deflation is repeated. 

The unitary DC method can be used in two ways. One approach is to 
divide the original e i gen prob 1 em , as well as subproblems so obtained, 
until only trivial e i gen prob 1 ems of orders two and one remain. fhese 
small e igenproblems are solved analytically. From the solutions of small 
e igenprob 1 ems , the solutions of e igenprob 1 ems of larger size are computed, 
and this step is repeated until the solution of the original eigenproblem 
has been determined. This approach is used in the numerical examples of 
Sect i on 5 . 

An alternative approach is to use the DC teclinique to generate just a 
few sube i genprob 1 ems , each of which can be solved independently by some 
other numerical scheme, sucli as tlie unitary GR method [Gr'l] , or the scheme 
in [AGRl] , in case the matrix is real orthogonal . 

We conclude this section with some bounds of the computational 
complexity of the unitary DC method. Assume that H G is given in 

Schur parametric form (1.1) with positive subdiagonal elements (Tj . Let n 
2 for some positive integer £, and subdivide the given eigenproblem until 
2 e i genprob 1 ems for 2x2 matrics are obtained. The latter e i genprob 1 ems 
are solved analytically. We assume that the number' of iterations required 
by the rootfinder for ^{0) can be bounded independently of n. 

Let first n independent pr^ocessors be available. Tire reduction of the 
original eigenproblem for H to R e i genprob 1 ems for 2x2 matrices can be 
carried out in t^ := 0(log22) finie steps. This computation only requires 
the determination of the Schur parameters for the unitary matrices of the 
smaller e i genprob 1 ems , see (2.3). Let the Schur’ par’ameters for all ~ 
unitary 2x2 matrices be known. The spectral resolution of all these 
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matrices can be computed in t 2 := 0(1) time steps. Assume that the 

partial spectral resolutions (2.21) of all unitary 2^ ^ x 2^ ^ matrices 

are known for some j G [2 , £] . In order to compute the partial spectral 

resolution of all 2^*^ unitary 2^ x 2^ matrices, we have to compute 2^ zeros 

£: 

of each of the 2 functions ^(^) , see (2.8) . Hence, a total number of n 
zeros have to be computed, and we use one processor to determine each one. 
Each function ^(^) has 2^ terms, and can therefore be evaluated in 0(2^) 
time steps for each value of 6. Hence, we can determine all eigenvalues 
of all 2^^ unitary 2^ x 2^ matrices in tj^ := 0(2^) time steps. For each 
eigenvalue we compute the first and last elements of the corresponding 
eigenvector from (2.21) . The first and last element of one eigenvector ^ 

can be determined by one processor in 0(2^) time steps. These 

computations have to be carried out for n eigenvectors by n processors and i 
therefore require t^^^ = 0(2^) time steps. Hence, the number of time steps 
required to determine the partial spectral resolution of H by n processors 
i s 






C 

+ ^2 + £ 
j=2 




€ 

+ E 

j=2 




0(n) 



(2.27) 



Now let n independent processors be available, and assume that the 
partial spectral resolutions of all 2^’^'^^ unitary 2^ ^ x 2^ ^ matrices are 
known for some j E [2 , €] . We have now n processors available for each 

£_j 

evaluation of each of the 2 functions ^(^). Each of these functions ^(^) 
has 2^ terms and can for each v'alue of 0 be evaluated in 0 (log 22 ^) time 
steps. Hence, we can compute all eigenvalues of all 2^^ unitary 2^ x 2^ 
matrices in tg^ := 0(log22*) time steps. The first and last elements of 
each eigenvector of each of the 2^^ unitary 2^ x 2* matrices can be 
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determined in := 0 ( 1 og 2 2^) time steps, by using n processors to 

compute each sum with 2^ < n terms. The initial determination of tlie ^ 

unitary 2x2 matrices and their spectral resolutions cannot be sped up 
essentially by using more than n processors, and requires 0 (t 2 ) -f 0 (t 2 ) 
time steps. Hence, the number of time steps required in order to compute 
the partial spectral resolution of H by n^ processors is 



0(ti) + 0 (t 2 ) + E + E = 0 (log 2 n) + E 0(j) = 0(log“n) . (2.28) 

j=2 j=2 j=2 



The time complexities (2 . 27) - (2 . 28) suggest tliat the unitary DC method 
presented could be attractive for use in real-time signal processing 
app 1 i cat ions. 



3. 



Some proper!, ies of tslie unitary DC met:hod 



We show some properties of the eigenvectors of II and zeros of 
Analogous results have previously been obtained by Dongarra and Sorensen 
[DS , Lemmas 4.2, 4.6 and 4.7] for the DC method for the eigenproblem for 
symmetric tridiagonal matrices. The following formulas are used in 
several of the proofs: 



^(A) = 1 + 21: iCj|2 

j=i ^ 



— 2 ^ 



AA) = -2 E IC+ 



(3.1) 



(3.2) 



<I-'(0) = ^'(A)A . 



(3.3) 



Lemma 3 . 1 Let A^/i E C, |A| = |/i| = 1 and A ^ /i , Assume that X^fi 0 {Aj}j|_p 
where A = diag[A 2 ,A 2 , . . . , Ap] . Let v^^ and v^ be defined by (2.20) . Then 



f r 1-1/2 



= |<A'(A)/(/i) I 



= |<I.'(^;^)<I.'(0^)| 



^(A) - 



(3.4) 



1/2 






\e^ \e^ 

e ~ e 



i^A *^/i 

where A = e ,/i = e 27 t . In particular, if A and /i are 

distinct eigenvalues of H, then <?^(A) = <^(/^) = 0, and therefore = 0. 

Proof. By (2.20), (2.6) and (2.7), 





’ Wl 




’ Ai 


V. = 








A 


W2 




I n-s 



(A - IA)-lz/6(A) , 



(3.5) 



and therefore 
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z 



V = 



(A - 1A)-1(A - 1/0 



A - f{X) 



-1 






= (6(A)6(/0)-' E 



Kjl 



i = l (Aj - A)(0 - /O 



(3.G) 



Now 



AA: 



(Aj-A) (Aj.-/i) (Aj-A)(Aj-/0 A-/x\Aj 



a_/_Al 

X-LlJ 



(3.7) 



Subst. i t<ut. i ng (3.7) int-o (3.6) yields 



v,"v, = (2«(A)H.))-‘ '4^' - 



2 I liilllS' 



j=l 



/x-A. 



and by (2.19), (3.1) and (3.3), 



= (.^'(a)^(/,)a/x)-i/2^ 



= ie 



1/2 . ^(^a) - 



i^A i^/x 

e - e 



This shows (3.4) . 



□ 



The denominat-or | A-// | in (3.4) suggests that it rna>^ be numerically 
diTTicult to obtain orthogonal eigenvectors when the associated 
eigenvalues are very close. The Tol lowing lemma sheds some light on this 
j situation, and shows that due to deflation the roots of ^’(^) are bounded 
away T rom each other. 

Lemma 3 . 2 Let A j = e ^ , 1 < j ^ be the eigenvalues of A, and let z = 

j [Cj] jl_i be defined by (2.7) . Let be an arbiti'ary but fixed positive 
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I 



constant and assume that the Aj are pairwise distinct and that 2 | (j | > 



tor all j. These conditions can be made valid by detlation. Assume that 
the Aj are sorted so that 0 < < ^2 • • • < < 27 t , and let 27 t + 

9^. Let 9 G [0,27t[ be a zero of ^(^), and let k be such that ^ • 

Then 



,2 

+ - 2^^k + l”^k) ^ ^"^k ^ ^^^^{fgC^k + l^^k) ’ ^} • (3*9) 



Proof. Introduce the index sets 



:= {j : ^ -f 27 t£ < ^ , for some £ E Z, 1 < j < n} , 



I 2 := { j : 9-7T < + 27t£ < 9^ for some £ E Z, 1 < j < n} 



Then D I 2 = 0 and 1 ^ U I 2 = {1,2,.. .,n}. Fu rt he r 



A + Aj i9 - (9.^ f < 9 , j E Ii, 

r-i-A = 



> 0 , 



j e i2- 



In particular, k E I 2 and, provided that k < n, kfl E I^. If k = n then 
1 E Ij. Moreover, 



cot 



( 2 ^) ^ ) ’ Vj G Ii 



cot 



(^) ^ . 



e - e..s 



(3.10) 



Vj G I. 



From =0, it follows that 



^ ~^l\ o i\ 



(3.11) 
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11 

By (3.10), (3.11) and z z = 1 we obtain, provided that k < n, 



-iCk + il^ cot( ^ ^2^^ ) ^ 
or, equivalently, 

iCk+ii^ ^ • (3. 



It k = n then we define Cp-fi •= Cn stud (3 . 1 2 ) - (3 . 1 3 ) remain valid. 
Now assume that 



^ ^ 2^^k + l “ ^k) • 



(3. 



We wish to determine a lower bound tor ^k + l”^* From tan ^ ^ tor 

0 < X < ^ it tollows that it 0 < ^k + l“^ — then 

-) ^ ^k + l - ^ ■ (3 

Substituting (3 . 1 4) - (3 . 1 5) into (3.13) yields 



ICk + ll^ + l “ ^k)) ^ ^k+1 “ ^ ’ 

and trom (^ (^^^ ^ ) ) > ^(^k + l“^k)' obtain 

I ^k-fl I ^ 5 (^k + 1 " ^k) ^ ^k + 1 “ ^ 

Finally, substituting Kk-fil ^ ^ into (3.16) yields (3.8). 

In order to show (3.9), we note that trom (3 . 1 0) - (3 . 1 1 ) and z^z = 
to 1 I ows that 

- co-fc( 2~^) - 

or , equ i valent 1 y , 

'tan(^ — 2“^) - tan^— ) , (3 



12 ) 



13) 



1-1) 



.15) 



16) 



1 i t 



21 



which corresponds to (3.13). We now assume that 



^k + i - ^ > J(^k+i - • (3.18) 

We would like to determine a lower bound for 0-6^. Similarly as in the 
derivation of (3.15) we obtain that if 0 < then 

tan( ^ ~ < e - . (3.19) 

From (3 . 17) - (3 . 19) and t^an (i (^j^_l_^-^)) > j^-^) , we obtain i 

^ > Kkl"i(^k + l - ^k) • (3.20)^ 

Finally, substituting I Ck I ^ ^ into (3.20) yields (3.9). □ 

Our final lemma shows that the computed eigenvectors are close to 
orthogonal if the zeros of ^(^) are evaluated with sufficient accuracy. 

Lemma 3.3 . Let A = diag[Aj,A 2 , . . . , , and let A, p, be computed 

approximations of the distinct roots A, fj, of <t> • Introduce the relative t 

errors , (3^ of Aj^-A and X^-p ^ respectively, i.e. 



jAk - A - (A^ - A)(l + a^) 

Ia^ - A = (Ak - ^) (1 + ^k) 



k = 1 , 2 , . . . , n 



( 3 . 21 ) 



Assume that for some constant 0<€<1, |<>kl — ^ and | /?k I — ^ for all k, 

and that |A| = \p\ = 1. Then 




|v;^''Cv^| < e(2 + e)(L4_i)2 



where C = d i ag [p j , » • • • » Pn] with 



22 



_ °k + + Qk^k ( «^(A)6(/i) A 

' (1 + «k) (1 + /?k)U(A)6(/i) 7 



(3.22) 



For Tj := 0 < < 27 t , we define 6(i/) := ( ^ 



1/2 



\l/2 



Proof. We first note that since A, /i are computed by determining zer^os of 
0 < 6 < 27t , the requirement |A| = \ fl \ = 1 is satisfied. Analogously 
to (3.6) we obtain 



A 



= {6{X)6{roy^ E 



1 

' E 




CM 

El. 




k=l 


(\ 


“ ('^k 


- /'O 


1 






Kkl'^ 


k=l 


(\ 


- A)(Ak 


- /0(i + “k)(l + ^k) 



where the last equality follows from (3.21) . Now 



0 = = ( 6 (A) 6 (, 0 )-' t 



Kkl" 



k=l (^k “ ''') (-^k “ /O 



and (2.19) imply that. 

Kkl" 



n 

E 



k=l ('^'k “ (-^k ■ f) 



= 0 . 



Theref ore 






n 

- E 



Kkl 



k=i (^k-^)(^k"/0 



= (sC^)s(fOr^ E 



ICkI 



ktl (Ak-A) (A^-,0Ul+5k) (l+/?k) 



(, 



(3.23) 
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which shows that ~ '^i'th C defined by (3.22) . 

(3.3), (3.2) and (3.21) it follows that 

^ ^'(A)a 
\S(X)J 4>'CX)X 





Kkl" 



(A-A,)" 



(-Ar)A VA n 

(A-A,)2 (i+a,)2>' 



F rom (2.19) , 



(3.24) 



F rom 



(-AkA)/(A-Aj 

VA 

(1 + «k)" ■ 



^ > 0 it follows that A/ ( A ( 1 +Qj^) > 0 and therefore 

(1 + kkl)'^ > (1 + O'^ • (3.25) 



Substituting (3.25) into (3.24) yields 6(A)/6(A) < 1 + e, and similarly one 

can show that <5(A^)/^(/i) <1+6. Hence, by (3.22), 

iPki < ^ (1 + o" = ^(2 + o(r4-7)' • (3-26) 

F i nal 1 y , 

- II^aIIs IICII2 \Wfih = IICII2 = l^kl ’ 

and the desired bound now follows from (3.26) . □ 
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4 . The ort:hogonal e i genprob 1 ein 



The computat: i odclI work required for The real orThogonal eigenprohlem 
is smaller Than Tor The uniTary one. 4'h i s secTion discusses These 
diTTerences, and considers some deTails of our i rnp 1 emenTaT i on of a DC 
scheme for The real orThogonal e i genprob 1 ein . Our compuTer program is for 
The case when 11 € wiTh n = 2^ , where £ is a posiTive iiiTeger, and we 

assume in This secTion ThaT n is of This form. Many of our commeriTs are 
valid for more general values of n, also. 



We f i rsT noTe ThaT The subdivision of The e i ge n pro b 1 ein for H i nTo 
I smaller e i genprob 1 ems , as described by (1.3) -(1.7), does not require any 
I compuTaT i onal work. Subdivision yields The block-diagonal maTrix 



H G1G3G5 . . . Gn.5Gn.3Gn.1Gn , 



(4.1) 



and we obTain simple formulas for The eigenpairs of each 2x2 block on The 
diagonal as follows. LeT 



G : = 



-7 

cr 



a 

7 



>2X2 



-1 <7< 1, (7>0, + (j^ = 1 



(4.2) 



Since G is real, symmeTric, orThogonal and has disTincT eigenv'alues 
{Ai,A 2 }, we have = 1 and A 2 = -1. I^eT eigenv'ector of 

uniT length. Then we can choose 



'^1 = 
.^2 = 



(4.3) 



Bind from a = ( 1 - 7^)^^^^ it follows ThaT 



_ ^)l/2 

^2 = (1 - 



(4.4) 
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Cance 1 1 at, i on of significant digits is avoided by using (4.3) if 7 > 0 and 
(4.4) otherwise. An eigenvector associated with A 2 = -1 is given by X 2 : = 

If 7 n = -1 then Gp = In 5 and the eigenpairs of G^.^Gn are those of Gp,_j . 
If 7 y^ = 1 we need to determine the eigenpairs of 



G' : = 



-7n-l 



'n-1 



^n-1 'yu-1 



We find that the eigenvalue = - 7 n-i + ^^n -1 has an associated 

eigenvector := ^ , and the eigenvalue A 2 = “Tn-i “ i^n -1 has an 

associated eigenvector X 2 := [ 2 '^^^ , ^ . 

Note that since the eigenvalues of G given by (4.2) are A = ±1, 
independent of -1 < 7 < 1 , deflation takes place numerous times during the 
computations . 



We turn to the computation of the Householder transformation (1.4) . 

In oder to avoid cancellation of significant digits, we compute {ws,Ws^i} 
given by (1.6) as follows. If 75 > 0, then we use (1.6a) and replace 
( 1 . 6 b) by 

Ws + I := - • ( 1 . 6 b') 

In case 75 < 0, we use (1.6b) and replace ( 1 . 6 a) by 

Ws := (Ts 2-1/2 , ( 1 . 6 a') 



Due to H having real-valued 
of H occur in complex conjugate 
0 < 6 < n have to be computed , 
n « /^i ~ ^ \ 



elements, the eigenvalues and eigenvectors 
pairs. Therefore only zeros of ^(^) for 
Moreove r , 
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can be simplified. Assume that 0 < < tt for some k < n and let 



= 27t - 0 ^ . Then | Ck I == I Ck-f i I ’ and we obtain 



Kkl^cot(^5^) + Kk + 1 = iCk - cot (5^)) 



-0s 



^O^-Os 



rO^+Os 



= ICki" 



2 sin 6 



cos $ 1 , I I 



sin 6 






(4.5) 



We use the right hand side of (4.5) in the computations, 
we need to evaluate cot(-g) as well as — tan 

The contribution from (4.5) to ^^(^) is 



2|Cu 



2 _cL(. 



sin 0 



d^?\cos 6 - cos 0, 



) = 2KJ 



2 1 - COS 0 cos 0^ 

(cos 6 - cos ^k)^ 



If = 0 then 



(4.6) 



The stable evaluation of the right hand side of (4.6) can be accomplished 
as described in Table 4.1. 



Cond i t i ons 



Eval uate 



CC|^ <0 



ccu > 0 and 



s^s- 



> 0 



cci, > 0 and g— < 0 



1 / -2 -2n 

) 



2s^s- 






( 2 *-.)= 



2s_^s- 






V2s^sJ 



Table 4.1: Stable evaluation of ( 1 - ccj^ ) / ( c - c^ ) ^ . where 

( 0 0 V 

— ) ■ 

s. sin(-^j. 



The interlacing of the zeros of (^^(A) w4th the {Aj^}|J_j implies that it 
easily can be determined w'hether 0 = 0 or 0 = n are zei^os of ^’(^) • Let 
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t.he 6 ^^ 1 < k < n, be ordered so 'that, 0 < ^ ^2 — • • * £ ^ ^ < ^P + 1 

< ... < < 2 ;r . Since t,he = exp(i^|^) appear in complex conjugate 

pairs, we obtain 

re^ > 0 => <i>( 0 ) = 0 , 

i^p < 7T => <I>(7r) = 0 . 



Finally, we consider the computation of eigenvectors defined 

by( 2 . 20 ). Let 



11 

1 — 1 




1 1 

CO 


J 


e c* , 


= : Cwf> 




‘ • 5 '^n-sj » 


wP) 

J 


e c"-® , 



Aj =: diag[exp(i^j^^) ,exp(i02^^) , . . . ,exp(i4^^)] 

(2) (2) (2) 

Aj =: diag[exp(i 0 j ),exp(i ^2 )>•••» exp ( i^n-s) ] 

=: [cS^4^...,ci¥ , 

W a, -• 

W2 e^Wg^j^ — . ,(^2 j • • • 9sn-sJ > 



0 < < 2 n , 

0 < < 27 t , 



and A =: exp(i^), 0 < ^ < 27 t . Then 

Wi(I-Ai^A)-lWi^esa,s = E (1 - exp ( i ) )'^ 

j=l J J J 

= E (1 - exp(i(0-^p4))'^Cp^w|^^ 

( 1 ) J J J 

^e{ 0 ,:r} 

+ E [(l-exp(i(0-0|^4))''Cj^^^^ + (l-exp(i(^+^j^4))-lrf^y^^] 

^1\ J JJ J JJ 

0 <^j 
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“2 (1 + i cot (|) + 2 ^ ( 1 - i tan (|) 



(9^*^=0 
J 






sin 



(1) 



^ V r^n rA^"> x 
+ E [(f^eCCj Wj ) + 

2 sin(-L_)si„(xL_) 






....«r»cr>>] 



(^1.7) 



qA'>-0 

- i sin e Y. ^sin^-j-2 — ^sin^^-^ — 



We may assume that close eigenvalues have been eliminated -from and 

(1) (2) 

by deflation, and that tlieretore tlie 0> and 6> are distinct. Hence, the 
sums over = 0 and ^ = tt contain at most one term each. 



Analogously to (4.7) we obtain 



V\^2(A2-lA)"^W2^eiu;s+i = E (^>^P ( ^ -^>^P ( i <^) 

j=l J J 



n-s 



o(2). 



.(2) (2) 



2 E (1+i cot(|) + i ^Y (-l + i tan (|) 



0^^^(2) (2) 






,( 2 ) 



,( 2 ) 



+ i t ^ 

J SI 



2 ) 2 ) 



o(2) 



/-'j A . ^ 

"( 2 ) ®^"( 2 ) 



n (2) 

(2) ) "j > 



Je-'^ Y sin 0j‘)(sin(— J-)sin(— ^))-l lm(Cj^tj^^) 



0<^j^^<7r 



+ ^ sin $ 



Y cos ^|2^(sin(^L)sin(^!-))-' 



, 0 + 0 . 



( 2 ) 



0 - 0 . 



( 2 ) 



(2) (2)^ 



0<f'<, 



(4.8-) 
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The evaluation oT 6(A) defined by (2.19) can also be simplified. We have 



(1) I 2 



^(A) = E 



!<■ 'I 



n-s 

+ E 






1 1/2 



i-l |A-exp(i6'j^^)|^ j-i I A-exp( i6lj^^) I 



where , e . g . , 



s 

E 



K 



( 1 ) 



i E lcPV/sin"(|) + J E |CJ'^lVcos2(|) 

4 J 2 4 J 



,( 1 ) 






j-1 I A-exp(i0j ’) I 






= TT 



+ 3 E Kj 
( 1 ) •' 
0<6lj ’<n 



( 1 ) 



6 0 ^ 0 ^^^ 



(4.9) 



The s i mp 1 i f i cal: i ons of “this section for the orthogonal eigenproblem have 
been implemented in a Pascal program. Several other mathematically 
equivalent forms of (4.7)-(4.9) could also be used. We have tried to find 
formulas that avoid unnecessary loss of significant digits. 
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5. 



Numer i cal examp 1 es 



We report results of some computed examples with an experimental 
program for the orthogonal e i genprob 1 em . Tiie program is written in Turbo 
Pascal 4.0 and was run on an IBM AT computer with unit roundoff 

u = 2*^^ ^ 2*10'^^. Our code implements the formulas of Section 4. 

Generally very accurate answers are obtained. Lemma 3.2 indicates, 
however, that a zero 9 of may be very close to a singular point 6 ^ of 

and by Lemma 3.3 the difference 9 - 9 ^ has to be computed to high 
re 1 at i ve accuracy in order to yield nearly orthogonal eigenvectors. 

Example 5.2 below shows that, indeed, 9 - 9 ^ can be extremely tiny and that 
loss of accuracy in both eigenvectors and eigenvalues may result. This 
loss of accuracy could be reduced, e.g. , by representing 9 and 9 ^ in liigher 
precision aritlimetic. 

’ In this section A G denotes the diagonal matrix witTi the computed 

eigenvalues of H G as entries, and W G is the matrix with the 

computed eigenvectors. We evaluate the residual errors ||HW-WA||oo and 

I I ||W*^W- 1 Iloo 9 where || ||qo denotes the uniform matrix norm. 

1 

I Example 5.1. This example discusses the application of tiie unitary and 
jlorthogonal e i gen prob 1 ems to the construction of Gauss-Szego quadrature 
jjrules. Consider the inner product on tiie unit circle 



<T,g> = [ f(A) g(A) da(A) 

J|A| = 1 



(5.1) 






vith a positive measure da (A) . Let 0 < k < n, be monic orthogonal 

polynomials with respect to (5.1). They satisfy a recuiM^ence relation 



1 
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r V’o(A) := 1 (5.1a) 

I := + TkV’k-i(A), 1 < k < n, (5.2b) 



for some parameters 7 (^ 6 C such 'that | Tk I ^ for 1 < k < n. Here 

^V’k-i ( 1/A) is the ’’reversed polynomial.” Let 7 n € C be an 
arbitrary complex number ot unit magnitude, and detine V’n t>y (5.2b) with 
k := n. Writing the recursions (5.2) (tor 1 < k < n) in matrix torm^ 

yields the unitary matrix 

H = G 1 G 2 . . .G^_iGn, (5.3) 



whose eigenvalues are the zeros ot t/’n • Here is detined by 7 (^ 

according to (1.2) tor 1 < k < n. Hence, the parameters { 7 ^}^ are the 

Schur parmaeters tor H. Let H = WAW be a spectral resolution, and detine 

T 2 

the weights := |e^Wej^| tor 1 < k < n. Then 

[ f(A)da(A) = t Pk^(\) + ^n(f) 

J|A|=1 k=l 



is a Gauss-Szego quadrature rule with respect to the measure da(A), 
because the error ^n(^) vanishes when t is any trigonometric polynomial 
degree less than n. See [Gr2] tor details. The computed examples 
illustrate the case when all Schur parameters 7 ^ are real valued and H 
theretore is real orthogonal . 

A particularly simple example is 7 ^^ := 0, 1 < k < n, and 7 n := ~1. 

Then V’k(A) = A^ , 0 < k < n, and V’n(^) = A" - 1 , and theretore 

= exp (27 t i (k - 1 ) /n) , 



1 < k < n . 



(5 



Pk 



1 /n . 



of 



4 ) 
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These Schur parameters have been used for Table 5.1. In the tah>le “# 
defl. close e.v.” stands for number of deflations due to close 
eigenvalues, and defl. small I I ” short for number of deflations 

due to components of z of small magnitude. "fwo eigenvalues are 
considered close if (2.26) is satisfied for ^2 •= 1*10’^, and 1 Ck I 

regarded small if (2.24) is valid for := 1-10"^. These values of and 
^2 are used in all computed examples of this section. 



n 


llHW-WAlloo 


llW^W-Illoo 


# defl. close e.v. 


# defl . 


smal ] 


4 


4 .6-10'^^ 


1 .510‘“ 


2 




0 


8 


6.4-10'^^ 


3.010'^^ 


8 




0 


16 


1 ,4-10'“ 


5.4-10'“ 


24 




0 


32 


2.9-10'^^ 


1 .810‘^° 


64 




0 


64 


3.910‘“ 


3.4-10'^° 


160 




0 



Table 5.1: : =z 0 , 1 < k < n ; 7n -1 

For 7 j^ :=0, 1 <k<n, and 7 p := 1, we obtain the polynomials = 

, 0 < k < n, and ^n(A) : = A^ -f 1. Hence, the eigenvalues are A^^ = 

exp ( i 7T (2k- 1 ) /n ) , 1 < k < n, and the Gauss-Szego weights are the same as 

in (5.4). Table 5.2 shows computations for the present Schur parameters, 
and differs from Table 5.1 mainly in that fewer deflations take place. 



n 


||HW-WAlloo 


l|W^W-l||oo 


# defl. close e.v. 


# defl . 


smal 1 


4 


7.8-10‘^^ 


1 .6-10‘“ 


0 




0 


8 


1 .710'“ 


4 .2-10’“ 


2 




0 


16 


3. 110'^^ 


1 . 610'^° 


10 




0 


32 


4 . 1 -10'“ 


3 . 5-10'^° 


34 




0 


64 


5 .6-10'^^ 


7.510'^° 


98 




0 






Table 5.2: 


Ik -=0, 1 <k<n; -)n 


:= 1 





I 
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In Tables 5. 3-5. 4 we have chosen 7 |^ := 0.8, 1 < k < n. This makes t,he 
A,^ gather in the left half plane. For the examples of Table 5.3 we have 



max Re 


< "i- 


For the examp 


les of Table 5. 


4 we 


obfai 


n max Re A|^ 


n 

4 


||HW-WA||oo 

2.7-10'^^ 


l|W^W-I||oo 
1 .610'^^ 


# def 1 . cl ose 
2 


e . V . 


# 


def 1 . smal 1 
0 


8 


5.5-10'^^ 


1 .810'^° 


8 






0 


16 


5 .2.10'“ 


3.2-10'^° 


24 






0 


32 


3. 2- 10'® 


9.3-10'® 


63 






1 


64 


3.2-10"® 


1 .6-10‘^ 


157 






3 




Tab 1 e 5.3: 7 ;^ 


VI 

00 

o 

II 


< n ; 


7n : = 


-1 


n 

4 


I|HW-WA||oo 
4 .8-10"“ 


l|W^W-I||oo 
1 .7-10‘^° 


# def 1 . c lose 
0 


e . V . 


# 


def 1 . smal 1 
0 


8 


9 .410'^^ 


5.5-10'^° 


2 






0 


16 


4 .8-10'^° 


6 .6-10'® 


10 






0 


32 


6.310‘^° 


2.310'® 


34 






0 


64 


4.2-10'® 


1 .9-10'^ 


97 






1 



Table 5.4: 7|^:=0.8,l<k<n;7n’=l 

In the last computed quadrature rules of this example we let the 7 |^ , 1 

< k < n, be uniformly distributed in the open interval and let 7n 

be -1 or 1 with probability ^ each. The 7 |^ are determined with the random 
number generator of Pascal. Table 5.5 shows the result of 30 
e i genprob 1 ems so generated. The maximum, average and minimum in Table 5.5 
are over all 30 e i genprob 1 ems . 
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IIHW-WAIloo 


l|W^W-l||oo 


# defl. close e.v. 


# defl. smal 1 | 


max 


7.2-10'^ 


2.5-10‘^ 


30 


0 


average 


5 . 8 - 10 '® 


2 .3-10‘^ 


26.5 


0 


min 


2.9-10'® 


1 .5-10'® 


22 


0 



Table 5.5: UniTormly distributed 7 ^^ 1 ^ ^ uniformly 

distributed 7 p E {-1,1}. Max, average and min are over 
30 e i gen prob 1 ems with n := 32 



I 



The numerical experiments of Table 5.5 indicate that for many choices 
of Schur parameters 7 ^^ , the magnitudes | Ck I slvg not sufficiently small to 
give rise to frequent deflations. This behavior' has also been observed in 
many other computed experiments. In contrast, massive deflation in DC 



I 



methods for symmetric tridiagonal matrices often is caused by small 

i components of the vector correspnding to z = [Ck] ^ 

( 

j Example 5.2. This example suggests that it might not be possible to 
I increase the small lower bound foi^ min | | of Lemma 3.2 significantly. 

i| The Schur parameters for Table 5.6 are obtained by reversing tire sign of 

I 

J the 7 |^ , 1 < k < n, of Table 5.4. 



1 n 


min |6'k| 
1 <k<n 


IIHW-WAIloo 


l|W^W-I||oo 


# defl . 
c 1 ose e.v. 


# defl 
sma 11 1 C 


4 


6.6-10'^ 


6 .910'“ 


3. 110'^° 


0 


0 


8 


00 

0 

1 


7.210'^° 


2 .5-10‘® 


2 


0 


16 


0* 


1 .2-10’^ 


3. I IO'® 


10 


0 


32 


0* 


7.210'^ 


1 .910‘^ 


34 


2 


1 64 


0* 


7.2-10'^ 


2. 6- 10'® 


97 


5 



Table 5.6: 7 ^^ := -0.8, 1 < k < n; 7 n := 1 . *The matrix has 

numerically the eigenvalue A = 1 o f multiplicity' two. 



35 



Because = 0.6 >0, 1 < k < n, t:he matrix H has distinct eigenvalues 

mathematically. Numerically two eigenvalues are so close that they are 
not distinguished with our present choice of €2 ~ A smaller value 

of € 2 ? such as C 2 = 1*10’^, gave in some numerical experiments larger 
residual errors |lHW-WA||oo or ||W^W-I||oo« ^ 
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